setwd("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/")
library(ChIPpeakAnno)
library(meta)
fn=list.files(pattern=".dat")
fn=fn[-grep("report", fn)]

case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
case_t=data.frame(chr=case$chr,start=as.numeric(case[,3]),end=as.numeric(case[,3])+1)
bed_case=toGRanges(case_t,format="BED",header=TRUE)

col_names=c("TWAS.db","con.file.name","case_in_region","case_not","con_in_region","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
for(i in 1:length(fn)){
file=read.table(fn[i],head=T,sep="\t")
file$P0=format(as.numeric(file$P0),scientific = F)
file$P1=format(as.numeric(file$P1),scientific = F)
file$TWAS.P=as.numeric(file$TWAS.P)
file_sig=file[file$TWAS.P<0.01,]
file_sig=file_sig[!is.na(file_sig$PANEL),]
file_sig$unitID=paste0("chr",file_sig$CHR,":",file_sig$P0,":",file_sig$P1)
file_sig=file_sig[!duplicated(file_sig$unitID),]

file_t=data.frame(chr=paste0("chr",file_sig$CHR),start=format(as.numeric(file_sig$P0),scientific = F),end=format(as.numeric(file_sig$P1),scientific = F))
bed_file=toGRanges(file_t,format="BED",header=TRUE)

ol1=findOverlapsOfPeaks(bed_case,bed_file)
tcase=as.data.frame(ol1$peaklist$`bed_case///bed_file`)
for(j in 1:4){
con.file=c("nobias_AShM.txt","nobias_AShM_total_reads10.txt","nobias_AShM_p.5.txt","nobias_AShM_total_reads10_p.5.txt")
con=read.table(con.file[j],head=F,sep="\t")
con_t=data.frame(chr=con[,1],start=as.numeric(con[,2]),end=as.numeric(con[,3])+1)
bed_con=toGRanges(con_t,format="BED",header=TRUE)

ol2=findOverlapsOfPeaks(bed_con,bed_file)
tcon=as.data.frame(ol2$peaklist$`bed_con///bed_file`)

result_tmp=data.frame(matrix(NA,1,ncol=8))
names(result_tmp)=col_names

a=dim(tcase)[1]
b=727-a
c=dim(tcon)[1]
d=dim(con_t)[1]-c
result_tmp[,1]=fn[i]
result_tmp[,2]=con.file[j]
result_tmp[,3]=a
result_tmp[,4]=b
result_tmp[,5]=c
result_tmp[,6]=d
result_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
result_tmp[,8]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,result_tmp)
}
}
result=result[-1,]
write.csv(result,"TWAS_enrichment_0.01.csv",quote = F,row.names = F)

 result1=result[result$con.file.name==con.file[1],]
 result1$case_not=result1$case_in_region+result1$case_not
 result1$con_not=result1$con_in_region+result1$con_not

metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=result1,sm="OR",studlab = TWAS.db) 
forest(metaor3)

######文章用图
setwd("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/")
result1=read.table("TWAS_enrichment_精神类.txt",head=T,sep="\t")
result1$case_not=result1$case_in_region+result1$case_not
 result1$con_not=result1$con_in_region+result1$con_not
metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=result1,sm="OR",studlab = statue) 
forest(metaor3)

meta_data=as.data.frame(metaor3)
meta_data=meta_data[order(meta_data$TE),]
meta_data$num=as.numeric(c(2,3,1,4,5))
library(ggrepel)
meta_data1=data.frame(statue=meta_data$studlab,TE=meta_data$TE,num=meta_data$num)
 meta_data2=meta_data1
 meta_data2$TE=0
 meta_data_line=rbind(meta_data1,meta_data2)
ggplot(meta_data,aes(TE,num))+geom_point(aes(size=event.e,color=pval))+geom_text_repel(aes(label = event.e))+
  scale_color_gradient(low = "red", high = "green")+geom_line(data = meta_data_line,aes(x=TE,y=num,group=statue))+
  labs(x="log(OR)",y="status")+scale_y_continuous(breaks = meta_data$num,labels = meta_data$studlab)+
  theme_bw()+theme(panel.grid.minor = element_blank())